Chip/RNA relation analysis

Setup

source("scripts/utils.R")
source("scripts/settings.R")

load("data/processed/rna_dynamics.RData")
load("data/processed/chip_dynamics.RData")
load(paste0("data/processed/h3k4me3_data_sperm_", sperm_data, ".RData"))

Balloon plots

plot_balloonplot(rna_dyns[names(rna_dyns) != "ND"], chip_dyns[names(chip_dyns) != "Absent"], limit=50, x_caption="RNA expression", y_caption="H3K4me3")

ggsave('output/rna_chip/rna_chip_balloon_plot.pdf', width=5, height=3)

plot_balloonplot(rna_dyns_timing, chip_dyns, limit=100, x_caption="RNA expression", y_caption="H3K4me3", remove_dyns=c("NonZygotic", "Absent"))

ggsave('output/rna_chip/rna_chip_timing_balloon_plot.pdf', width=5, height=3)

rna_small <- rna_dyns[!names(rna_dyns) %in% c("ND")]
chip_small <- chip_dyns[!names(chip_dyns) %in% c("Absent")]
genes_of_interest <- unique(c(unlist(rna_small), unlist(chip_small)))
cols = lapply(c(rna_small, chip_small), function(dyn) genes_of_interest %in% dyn)
df <- data.frame(cols)
metadata <- data.frame(set = c(names(rna_small), names(chip_small)), type = c(rep("rna", length(rna_small)), rep("chip", length(chip_small))))
upset(
    df,
    colnames(df),
    name="dynamic",
    width_ratio=0.15,
    stripes=upset_stripes(
        mapping=aes(color=type),
        colors=c(
            'rna'='green',
            'chip'='orange'
        ),
        data=metadata
    )
)

Peak line plots

for(stage in stages){
  plot_peak(data=chip_norm_mats, stage=stage, dyns=rna_dyns, cols=cols_rna, remove_NA=TRUE, onlypeaks=FALSE, label = "H3K4me3")
  ggsave(file=paste0("output/rna_chip/", stage, "_rna_peak_plot.pdf"), width = 5, height = 4)
}

for(stage in stages){
  plot_peak(data=chip_norm_mats, stage=stage, dyns=rna_dyns_timing, cols=cols_rna_timing, remove_NA=TRUE, onlypeaks=FALSE, label = "H3K4me3")
  ggsave(file=paste0("output/rna_chip/", stage, "_rna_timing_peak_plot.pdf"), width = 5, height = 3)
}

Euler / Venns of overlaps

genesAbsolute <- lapply(chip_dyns, function(x){as.data.frame(lapply(rna_dyns, function(y){length(intersect(x, y))}))})
genesAbsolute <- do.call(rbind, genesAbsolute)

# all genes
plot(venn(c("Kept"=length(chip_dyns$Kept), "GZ"=length(rna_dyns$GZ), "Kept&GZ"=length(intersect(chip_dyns$Kept, rna_dyns$GZ)))),
     quantities = TRUE,fills=c(cols_chip["Kept"], cols_rna["GZ"], "#fa193b"))

plot(euler(c("Kept"=length(chip_dyns$Kept), "GZ"=length(rna_dyns$GZ), "Kept&GZ"=length(intersect(chip_dyns$Kept, rna_dyns$GZ)))),
     quantities = TRUE,fills=c(cols_chip["Kept"], cols_rna["GZ"], "#fa193b"))

# only genes where both chip & rna dynamics are known
euler_kept_gz <- function() {plot(euler(c("Kept"=sum(genesAbsolute["Kept",])-sum(genesAbsolute["Kept","GZ"]), "GZ"=sum(genesAbsolute[,"GZ"])-sum(genesAbsolute["Kept","GZ"]) , "Kept&GZ"=sum(genesAbsolute["Kept","GZ"]))),
     quantities = TRUE,fills=c(cols_chip["Kept"], cols_rna["GZ"], "#fa193b"))}
plot(venn(c("Kept"=sum(genesAbsolute["Kept",])-sum(genesAbsolute["Kept","GZ"]), "GZ"=sum(genesAbsolute[,"GZ"])-sum(genesAbsolute["Kept","GZ"]) , "Kept&GZ"=sum(genesAbsolute["Kept","GZ"]))),
     quantities = TRUE,fills=c(cols_chip["Kept"], cols_rna["GZ"], "#fa193b"))

euler_kept_gz()

pdf('output/rna_chip/kept_gz_overlap_euler.pdf', width=4, height=3)
euler_kept_gz()
dev.off()

## quartz_off_screen 
##                 2

Stacked bar of H3K4me3 ratios

cols_mod <- c(cols_chip[!names(cols_chip)%in%c("Recovered", "Lost", "Gained")], "Other"="darkgrey")
dyns_mod = chip_dyns[c("Kept", "Absent")]
dyns_mod[["Other"]] = unname(unlist(chip_dyns[c("Lost", "Gained")]))  #Recovered in other

genesAbsolute_mod = genesAbsolute[c("Kept", "Absent"),]
genesAbsolute_mod["Other", ] <- colSums(genesAbsolute[c("Lost", "Gained"),]) #Recovered in other

ggdf1 <- melt(cbind(genesAbsolute_mod, rownames(genesAbsolute_mod)))

## Using rownames(genesAbsolute_mod) as id variables

names(ggdf1)<- c("H3K4me3", "RNA", "count")
ggdf1$H3K4me3 <- factor(ggdf1$H3K4me3, levels=c("Kept", "Other", "Absent"))
ggplot(ggdf1, aes(fill=H3K4me3, y=count, x=RNA)) + 
  geom_bar(position="fill", stat="identity") + 
  theme_bw() + 
  ylab(NULL) + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), 
        axis.title.x = element_blank(), 
        legend.title=element_blank(), 
        legend.position="right") + 
  scale_fill_manual(values = cols_mod) + 
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  theme_pubr()

ggsave('output/rna_chip/ratio_h3k4me3_rna_dynamics_stacked.pdf', width=5, height=3)

Violin plots of H3K4me3 for RNA groups

for(stage in stages){
  g <- plot_stage_dynamics(chip_dfs, dyns=rna_dyns, cols=cols_rna, timepoint=stage, type="promoter.mean", plot = "violin")
  grid.arrange(g)
  ggsave(file=paste0("output/rna_chip/", stage, "_rna_dynamics_violin_plot.pdf"), g, width = 5, height = 4)
}

## [[1]]
## [1] "GS" "GZ"
## 
## [[2]]
## [1] "ZS" "GS"
## 
## [[3]]
## [1] "ND" "ZS"
## 
## [1] 7.960691 7.076170 6.191648

## [[1]]
## [1] "GS" "GZ"
## 
## [[2]]
## [1] "ZS" "GS"
## 
## [[3]]
## [1] "ND" "ZS"
## 
## [1] 7.026371 6.245663 5.464955

## [[1]]
## [1] "ZS" "GZ"
## 
## [[2]]
## [1] "GS" "ZS"
## 
## [[3]]
## [1] "ND" "GS"
## 
## [1] 5.091035 4.525364 3.959694

## [[1]]
## [1] "ZS" "GZ"
## 
## [[2]]
## [1] "GS" "ZS"
## 
## [[3]]
## [1] "ND" "GS"
## 
## [1] 6.624786 5.888699 5.152611

for(stage in stages){
  g <- plot_stage_dynamics(chip_dfs, dyns=rna_dyns, cols=cols_rna, timepoint=stage, type="peak.width", plot = "violin")
  grid.arrange(g)
  ggsave(file=paste0("output/rna_chip/", stage, "_peak_width_rna_dynamics_violin_plot.pdf"), g, width = 5, height = 4)

    g <- plot_stage_dynamics(chip_dfs, dyns=rna_dyns, cols=cols_rna, timepoint=stage, type="peak.width", plot = "ecdf")
  grid.arrange(g)
  ggsave(file=paste0("output/rna_chip/", stage, "_peak_width_rna_dynamics_ecdf_plot.pdf"), g, width = 5, height = 4)
}

## [[1]]
## [1] "GS" "GZ"
## 
## [[2]]
## [1] "ND" "GS"
## 
## [[3]]
## [1] "ZS" "ND"
## 
## [1] 5130 4560 3990

## [[1]]
## [1] "GS" "GZ"
## 
## [[2]]
## [1] "ND" "GS"
## 
## [[3]]
## [1] "ZS" "ND"
## 
## [1] 10955.7  9738.4  8521.1

## [[1]]
## [1] "GS" "GZ"
## 
## [[2]]
## [1] "ND" "GS"
## 
## [[3]]
## [1] "ZS" "ND"
## 
## [1] 7626.6 6779.2 5931.8

## [[1]]
## [1] "ZS" "GZ"
## 
## [[2]]
## [1] "GS" "ZS"
## 
## [[3]]
## [1] "ND" "GS"
## 
## [1] 8020.8 7129.6 6238.4

for(stage in stages){
  g <- plot_stage_dynamics(chip_dfs, dyns=rna_dyns_timing, cols=cols_rna_timing, timepoint=stage, type="promoter.mean", plot = "violin")
  grid.arrange(g)
  ggsave(file=paste0("output/rna_chip/", stage, "_rna_timing_violin_plot.pdf"), g, width = 5, height = 4)
}

## [[1]]
## [1] "Expressed@8hpf" "Expressed@6hpf"
## 
## [[2]]
## [1] "Expressed@7hpf" "Expressed@8hpf"
## 
## [[3]]
## [1] "Expressed@9hpf" "Expressed@7hpf"
## 
## [[4]]
## [1] "NonZygotic"     "Expressed@9hpf"
## 
## [1] 7.960691 7.076170 6.191648 5.307127

## [[1]]
## [1] "Expressed@8hpf" "Expressed@6hpf"
## 
## [[2]]
## [1] "Expressed@7hpf" "Expressed@8hpf"
## 
## [[3]]
## [1] "Expressed@9hpf" "Expressed@7hpf"
## 
## [[4]]
## [1] "NonZygotic"     "Expressed@9hpf"
## 
## [1] 7.026371 6.245663 5.464955 4.684247

## [[1]]
## [1] "Expressed@7hpf" "Expressed@6hpf"
## 
## [[2]]
## [1] "Expressed@8hpf" "Expressed@7hpf"
## 
## [[3]]
## [1] "Expressed@9hpf" "Expressed@8hpf"
## 
## [[4]]
## [1] "NonZygotic"     "Expressed@9hpf"
## 
## [1] 5.091035 4.525364 3.959694 3.394023

## [[1]]
## [1] "Expressed@8hpf" "Expressed@6hpf"
## 
## [[2]]
## [1] "Expressed@9hpf" "Expressed@8hpf"
## 
## [[3]]
## [1] "Expressed@7hpf" "Expressed@9hpf"
## 
## [[4]]
## [1] "NonZygotic"     "Expressed@7hpf"
## 
## [1] 6.624786 5.888699 5.152611 4.416524

for(stage in stages){
  g <- plot_stage_dynamics(chip_dfs, dyns=rna_dyns_timing, cols=cols_rna_timing, timepoint=stage, type="peak.width", plot="violin")
  grid.arrange(g)
  ggsave(file=paste0("output/rna_chip/", stage, "_peak_width_rna_timing_violin_plot.pdf"), g, width = 5, height = 4)

  g <- plot_stage_dynamics(chip_dfs, dyns=rna_dyns_timing, cols=cols_rna_timing, timepoint=stage, type="peak.width", plot="ecdf")
  grid.arrange(g)
  ggsave(file=paste0("output/rna_chip/", stage, "_peak_width_rna_timing_ecdf_plot.pdf"), g, width = 5, height = 4)
}

## [[1]]
## [1] "Expressed@7hpf" "Expressed@6hpf"
## 
## [[2]]
## [1] "Expressed@8hpf" "Expressed@7hpf"
## 
## [[3]]
## [1] "Expressed@9hpf" "Expressed@8hpf"
## 
## [[4]]
## [1] "NonZygotic"     "Expressed@9hpf"
## 
## [1] 5130 4560 3990 3420

## [[1]]
## [1] "Expressed@8hpf" "Expressed@6hpf"
## 
## [[2]]
## [1] "Expressed@7hpf" "Expressed@8hpf"
## 
## [[3]]
## [1] "Expressed@9hpf" "Expressed@7hpf"
## 
## [[4]]
## [1] "NonZygotic"     "Expressed@9hpf"
## 
## [1] 10955.7  9738.4  8521.1  7303.8

## [[1]]
## [1] "Expressed@7hpf" "Expressed@6hpf"
## 
## [[2]]
## [1] "Expressed@8hpf" "Expressed@7hpf"
## 
## [[3]]
## [1] "Expressed@9hpf" "Expressed@8hpf"
## 
## [[4]]
## [1] "NonZygotic"     "Expressed@9hpf"
## 
## [1] 7626.6 6779.2 5931.8 5084.4

## [[1]]
## [1] "Expressed@8hpf" "Expressed@9hpf"
## 
## [[2]]
## [1] "Expressed@7hpf" "Expressed@8hpf"
## 
## [[3]]
## [1] "Expressed@6hpf" "Expressed@7hpf"
## 
## [[4]]
## [1] "NonZygotic"     "Expressed@6hpf"
## 
## [1] 8020.8 7129.6 6238.4 5347.2